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Abstract 

The  effects  of  molecular  weights  (MAt  MB)  on  the  self-organized  segregation  of  immiscible 
constituents  ( A,B )  driven  by  pressure  bias  ( H  =  0.0- 1.0)  generated  by  geologic  processes  are 
examined  by  an  interacting  lattice  gas  Monte  Carlo  simulation.  Constituents  (A,B)9  released 
from  a  source  at  the  bottom  according  to  their  lattice  concentrations,  can  escape  the  lattice 
from  top  or  bottom.  The  longitudinal  steady-state  density  profiles  (AyB)  depend  on  their 
molecular  weight  and  bias  with  linear,  exponential,  and  non-monotonic  decays  with  the 
height.  The  transverse  density  profiles  show  signatures  of  partial  layering.  The  total  number 
(N)  of  constituents  remains  nearly  constant  in  steady  state  and  seems  to  decay  with  the  bias, 
N  oc  /r0-14  for  all  but  extreme  range  of  H.  At  MA  =  0.1,  MB> 0.8,  number  of  constituents 
( N  AyNB )  show  non-monotonic  dependence  on  H  with  NBPNA  but  a  maximum  in  NA  and 
minimum  in  NB. 
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1.  Introduction 

Flow  of  fluid  mixtures  from  a  reservoir  source  driven  by  the  fluid  pressure  bias 
resulting  from  compaction  of  sediments,  tectonic  processes,  geothermal  anomalies, 
etc.  is  a  common  natural  phenomena  in  geomarine  systems  [1—3].  Some  examples 
include  dewatering  of  fine-grained  sediments  due  to  overburden  pressure,  dissocia¬ 
tion  of  methane  and  hydrocarbon  below  the  ocean  floor  [4,5],  mud  volcanoes,  free- 
methane  and  water,  dispersion  of  complex  fluid  mixtures  consisting  of  hazardous 
waste,  etc.  Understanding  the  distribution  of  components  would  be  valuable, 
particularly  when  the  constituents  are  dissimilar  with  different  molecular  weights  as 
many  of  these  systems.  A  variety  of  patterns  [6]  have  been  observed  in  many  different 
mixtures  including  granular  systems  by  simple  computer  simulation  models  [7,8], 
e.g.,  pattern-forming  segregation  [9],  lane  formation  [10],  segregation  and  mixing 
[11-16].  Enormous  efforts  have  been  made  in  recent  years  to  understand  the 
structural  evolutions  in  particulate  mixtures  by  lattice  gas  [17,18],  molecular 
dynamics  [19,20],  and  Monte  Carlo  [21]  methods.  Distribution  of  dissimilar 
constituent  particles  leads  to  random  yet  well-defined  structural  patterns  [22,16]  as 
a  function  of  parameters  such  as  interaction  and  temperature. 

Why  the  distribution  of  constituents,  their  specific  mixtures,  and  underlying 
patterns  in  a  steady-state  flow  of  multi-component  system  are  important?  It  is  now 
well  known  among  the  geophysics,  geology,  oil  and  gas  exploration  and  research 
community  that  the  existence  of  methane  hydrate  below  the  ocean  floor 
nearly  throughout  the  continental  shelf  offers  abundance  of  energy  resources  that 
exceed  nearly  all  the  natural  and  fossil  fuel  together  [1,2].  The  kinetics  of 
hydrate  formation  [1],  i.e.,  the  chemistry  and  its  thermodynamic  stability  is 
well  understood  experimentally  and  theoretically  (computationally  in  particular). 
The  origin  of  methane  gas  below  the  ocean  floor  also  is  reasonably  well  understood. 
Reduction  of  organic  matter  by  bacteria  (generally  within  tens  of  meters 
of  the  seafloor)  generates  what  is  referred  to  as  biogenic  methane;  this  is  the 
most  widely  distributed  form  of  methane  in  the  marine  environment.  Leakage  of 
methane  from  deep  petroleum  reserves  through  faults  in  the  sediments  places 
thermogenic  methane  near  the  seafloor  especially  in  areas  with  significant  petroleum 
reserves  (e.g.,  the  Gulf  of  Mexico).  The  methane  gas  is  buoyant  and  flows  upward 
through  the  permeable  sediments  towards  the  water.  The  methane  is  under 
hydrostatic  pressure  and  at  a  temperature  determined  by  the  local  geothermal 
gradient;  when  the  free  methane  gas  (and  water)  reach  the  appropriate 
pressure-temperature  regime  methane  hydrate  can  form.  The  continuous  supply  of 
methane  gas,  formation  of  methane  hydrate  and  its  dissociation  are  on-going 
geothermal  processes  which  lead  to  their  steady-state  distribution  and  instabilities 
with  varying  time  and  spatial  scales.  Of  course,  the  enormous  magnitude  of 
scale,  diversity  of  emerging  phenomena  (dissociation  of  methane  gas  to  eruption 
of  mud  volcanoes  below  the  ocean  floor),  and  their  implications  and  impacts 
are  overwhelming  to  tackle  in  one  computer  simulation  modeling.  However, 
it  is  possible  to  examine  patterns  relevant  in  part  of  such  complex  system  and  learn 
the  basic  cause  and  effects  of  such  parameters  as  the  pressure  bias  and  molecular 
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weight  of  constituents  on  their  distribution  via  idealized  computer  simulation 
modeling. 

In  this  article,  we  consider  the  flow  of  immiscible  mixture  of  two  components 
(A,B)  of  molecular  weights  MA ,  MB  in  an  effective  medium  on  a  cubic  lattice  and 
examine  the  effects  of  molecular  weight  and  pressure  bias  on  the  steady-state  density 
distributions  by  an  interacting  lattice  gas  simulation.  In  particular,  we  address  how 
the  overall  steady-state  density  of  constituent  particles  and  their  profiles  depend  on 
the  pressure  bias  and  their  molecular  weights. 


2.  Model 

We  consider  the  mixture  of  two  components  A  and  B  (each  with  molecular  weight 
M a  and  Mb)  on  a  cubic  lattice  of  size  L3  with  L  =  30-200.  Particles  A  and  B  are 
randomly  distributed  initially  at  about  50%  of  the  lattice  sites  with  one  particle  at  a 
site  at  any  time.  A  nearest-neighbor  interaction  between  particles  (A,B)  and  empty 
(pore)  sites  ( O )  is  described  by 

£  =  E  E  ,  0) 

i  n 

where  index  i  runs  over  all  sites  occupied  by  particles  and  n  over  all  nearest-neighbor 
sites  of  i.  The  interaction  matrix  elements 

J(A,  A)  =  J(B ,  B)  =  - J(A ,  B)  =  - J(B ,  A)  =  -b  ,  (2a) 

J(AyO)^J(B}0)  =  - 1.  (2b) 

The  interaction  strength  e  =  1  describes  the  miscibility  gap.  The  gravitational  energy 
U  of  each  particle  with  mass  MAjb  at  the  height  2  is  given  by 

U  —  MAjBz  .  (3) 

A  source  of  particles  is  connected  to  the  bottom  plane  (z  —  1).  This  causes  a 
concentration  gradient  (see  below)  which  implicitly  exerts  a  driving  force  upward. 
Additionally,  we  consider  a  pressure  bias  (H)  caused  by  the  hydrostatic  pressure 
gradient  to  move  constituent  along  vertical  (±z)  direction  with  probabilities 

H+z  =  (1  +  #)/6,  H~z  =  (1  —  H)/6,  (4) 

Note  that  the  hydrostatic  pressure  bias  is  implemented  probabilistically  (see  below) 
and  acts  against  sedimentation  (Eq.  (3)). 

Each  randomly  selected  particle  (at  site  1)  attempts  to  select  a  neighboring  site  (/) 
along  ±x  and  ±y  with  equal  probability  (£)  and  along  z  (up)  and  -2  (down) 
directions  with  unequal  probabilities  (Eq.  (4))  for  non-zero  hydrostatic  bias  H.  The 
Metropolis  algorithm  is  used  to  move  randomly  selected  particles  (A  and  B)  to  their 
neighboring  empty  sites  with  probability  exp{-(AZs  +  AC/)}.  Note  that  the 
sedimentation  probability  is  coupled  with  the  gravitational  potential  energy 
(Eq.  (3))  of  constituents  via  Boltzmann  distribution.  When  a  particle  at  the  bottom 
plane  ( z  =  1)  moves  horizontally  (xy-plane)  or  upward  (to  z  =  2  plane),  a  new 
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particle  (A  or  B)  from  the  source  is  injected  onto  the  vacated  site  according  to  their 
current  lattice  concentrations.  Unlike  our  previous  studies  [15]  where  particles  A  and 
B  enter  the  lattice  with  equal  probability,  the  injection  of  particles  depends  on  their 
concentrations  (pA,pB)  here.  This  introduces  a  concentration  dependent  interaction 
between  the  source  and  the  system  which  lead  to  non-linear  responses  as  we  may  see 
below.  The  cubic  box  is  open  along  vertical  boundaries,  i.e.,  a  particle  can  drop  out 
when  it  attempts  to  move  down  from  the  bottom  or  escape  from  the  top  ( z  =  L). 
Periodic  boundary  conditions  are  implemented  along  the  transverse  ( x,y )  directions. 
An  attempt  to  move  each  particle  once  is  defined  as  one  Monte  Carlo  step 
(MCS)  time.  Almost  all  data  presented  here  are  generated  on  1003  sample.  Different 
sample  sizes  are,  however,  used  to  check  for  the  finite  size  effects.  No  severe  finite 
size  effect  is  observed  on  the  qualitative  behavior  unless  specified  otherwise. 
The  molecular  weight  ratio  a  =  MB/MA  is  varied  for  a  constant  MA  =  0.1  and 
a  =  1-10. 

The  number  and  densities  of  particles  and  their  distributions  change  as  the 
simulation  proceeds.  Injection  of  particles  from  the  source  at  the  bottom  leads  to  a 
net  flow  along  the  longitudinal  direction.  The  competing  driving  fields  due  to 
hydrostatic  pressure  bias  and  gravity  further  affect  the  flow  and  pattern.  A  steady 
state  is  reached  with  a  stable  density  profile  and  a  constant  rate  of  mass  flux  in 
asymptotic  time  limit.  The  simulation  is  repeated  for  a  number  of  independent 
samples  to  obtain  a  reliable  estimate  of  physical  quantities.  These  quantities  include, 
root  mean  square  (rms)  displacement  of  tracers  (particles)  and  that  of  their  center  of 
mass,  density  profiles  (longitudinal  and  transverse),  correlation  profiles  (i.e.,  average 
number  of  different  neighbors  of  each  particles),  flow  rate  and  flux  rate  density,  etc., 
as  a  function  of  bias  H  for  different  values  of  molecular  weight  ratios  a.  In  this 
article,  we  constrain  to  structural  response,  i.e.,  changes  in  density  profiles  of  A  and 
B  and  their  steady-state  concentrations. 


3.  Results  and  discussion 

Visualizations  (on  smaller  samples)  show  how  the  two  components  re-arrange  and 
segregate  to  a  steady-state  distribution  as  they  continue  to  flow.  Fig.  1  shows  the 
snap-shots  of  their  steady-state  distributions  for  different  bias  (H)  at  the  molecular 
weights  Mb  =  0.8, 1.0  for  MA  =0.1.  Sedimentation  and  segregation  are  clearly  seen 
at  lower  bias  (H  —  0.2, 0.3)  and  higher  a.  At  low  bias  (H  =  0.2),  sedimentation 
dominates  and  the  heavier  constituents  (B)  settle  down  more  toward  the  lower  part 
of  the  lattice  which  has  much  higher  density  than  the  upper  part  with  a  dilute 
concentration  of  lighter  constituents  (A).  Note  that  the  distribution  of  particles 
remains  the  same  while  both  constituents  continue  to  flow  from  bottom  to  top.  One 
may  interpret  the  high  density  at  the  bottom,  a  representative  of  dissociating  solid 
phase  while  the  dilute  region  toward  the  top  as  a  gas  phase.  Higher  bias  smears  the 
phase  separation  by  reducing  the  segregated  gaseous  domains  of  lower  molecular 
weight  component  (A)  by  interspersing  into  the  heavier  component.  Systems  remain 
segregated  with  regions  of  partial  layers  and  domains.  The  total  density  (number  of 
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Fig.  1 .  Immiscible  A  (light,  MA  =  0.1)  and  B  (dark,  MB  =  0.8  (left)  1.0  (right))  particles  with  pressure  bias 
E  —  0.2, 0.3, 1.0  at  time  steps  t  =  2500  with  e  =  1.  Injection  of  particles  is  proportional  to  their  lattice 
concentration. 


particles)  in  steady-state  depends  on  the  bias  and  molecular  weight  as  will  be  seen 
below. 

Fig.  2  shows  a  typical  mass  transfer  (bottom  to  top)  with  time  steps  as  a  function 
of  hydrostatic  pressure  bias  ( H  =  0.10-1.0)  for  a  =  8.  The  flux  rate  (the  slope  of  the 
mass  transfer  with  time  step)  becomes  constant  which  implies  that  our  system  has 
reached  the  steady  state,  i.e.,  the  number  of  particles  entering  the  system  from  the 
source  is  equal  to  their  number  leaving.  Accordingly,  the  total  number  of  particles  in 
the  sample  remains  constant  with  time  step  (see  the  inset  figures).  At  low  bias  (i.e., 
H  ~  0.0-0. 3),  the  flux  of  heavier  component  ( B )  is  lower  than  that  of  the  lighter 
constituents  (A).  At  high  bias  ( H  ->  1),  the  flux  of  B  component  becomes  much 
larger  than  that  of  A  resulting  in  a  larger  concentration  of  B  in  the  lattice.  This  is  due 
to  concentration-dependent  interaction  between  the  source  and  the  system-in-flux  of 
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Fig.  2.  Mass  transfer  versus  time  steps  for  various  bias  H  —  0.10-1.00.  Inset  figures  shows  the  number  of 
corresponding  particles  in  the  simulation  box.  Sample  size  1003,  Af^=0.1,  Mb  =  0.8  with  128 
independent  runs. 


particles  B  increases  (in-flux  of  A  decreases)  according  to  their  lattice  concentration. 
The  concentration-dependent  in-flux  from  the  source  introduces  a  long-range  (visco¬ 
elastic  type)  hydrodynamic  interaction  that  may  cause  volcanic  non-linear  flow 
response,  a  topic  of  out  future  analysis. 

As  mentioned  above  the  concentration  gradient  and  the  pressure  bias  are  driving 
the  constituents  against  gravitational  sedimentation.  In  steady  state,  the  density 
profiles  remain  constant  in  time  while  constituents  continue  to  flow.  Fig.  3  shows 
density  profiles  at  bias  if  =  0.1  and  0.4  for  different  values  of  molecular  weight 
ratios  a  =  1-4.  The  longitudinal  density  profiles  are  very  sensitive  to  a  at  low  bias 
( H  —  0.1)  where  gravity  has  profound  effect  leading  to  much  larger  densities  of  both 
constituents  ( A ,  B)  towards  the  bottom.  Still  there  is  a  net  mass  flow  from  bottom  to 
top.  The  transverse  (y)  density  profiles  show  oscillations,  a  signature  of  lane 
formation  or  partial  longitudinal  layering.  Note  that  the  complementary  spatial 
oscillation  in  transverse  density  profiles  of  A  and  B  constituents  implies  that  they  are 
segregated.  Significant  drop  in  the  longitudinal  density  profile  from  bottom  to  top 
also  shows  a  phase  separation  between  high  density  at  bottom  (mostly  B  component) 
and  low  density  A  (gas)  phases  (see  also  Fig.  1).  At  higher  bias  ( H  =  0.4),  such  a 
visible  layered  segregation  is  reduced  somewhat  as  the  A  constituents  interspersed 
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Fig.  3.  Longitudinal  and  transverse  (inset)  density  profiles  at  bias  H  —  0.1  and  0.4  for  M a  —  0.1  and 
Mb  =  0. 1-0.4.  Open  symbols  represents  A  and  filled  symbols  B  particles  densities.  Statistics  is  the  same  as 
Fig.  2. 


into  B.  The  linear  decay  region  of  the  longitudinal  density  profile  spans  with  a  larger 
density  of  higher  molecular  weight  component  (see  Fig.  3).  Thus,  the  steady-state 
segregation  of  immiscible  constituents  ( A ,  B)  and  phase  separation  of  their  high-  and 
low-density  phases  along  with  a  simultaneous  flow  show  dissociating  characteristics 
of  complex  evolutionary  systems,  i.e.,  the  formation  of  methane  hydrate  and  its 
dissociation  below  the  ocean  floor  [1-4]. 

Such  an  open  dynamic  system  in  which  the  number  of  constituents  are  not 
conserved,  exhibits  stable  density  profiles  and  distinct  phases  with  specific 
concentration  distributions  that  depend  on  the  hydrostatic  pressure  and  molecular 
weight  at  a  fixed  temperature  (T  —  1  in  unit  of  Boltzmann  constant  here).  Obviously 
the  total  number  of  constituents  (i.e.,  their  overall  densities)  depends  on  pressure  bias 
and  a.  Fig.  4  shows  the  dependence  of  the  number  of  particles  on  the  bias  H  for  a 
range  of  molecular  weight  ratios,  a  =  1—10.  We  immediately  note  a  symmetry  in  the 
distributions  of  these  constituents.  At  a  =  1,  the  distribution  of  A  and  B  coincides  at 
all  bias  with  nearly  a  constant  value  of  their  concentrations.  At  higher  values  of  a, 
the  number  of  B  constituents  (higher  molecular  weight)  become  larger  than  that  of  A 
with  a  highly  nonlinear  rate  of  change  at  extreme  values  of  H  cz  0.8-1. 0.  The 
variation  of  the  steady  state  (i.e.,  nearly  constant)  numbers  of  A  and  B  with  the  bias 
becomes  non-monotonic  at  high  molecular  weight  ratio  (a  =  8, 10).  While  the 
number  of  A  particles  shows  a  maximum,  the  number  of  B  particles  shows  a 
minimum  around  H  ~  0.83. 
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Fig.  4.  Number  of  A,B ,  and  their  total  numbers  in  steady-state  versus  bias  H  for  MA  =0.1  and 
Mb  -  0. 1-1.0  (upper  legend).  Top  inset  figure  shows  the  total  number  of  particles  versus  H  on  a  log-log 
scale.  Lower  inset  figure  shows  the  variation  of  average  transverse  (y)  density  with  Mb  at  H  —  0-1— 0.8 
(lower  left  legend).  Open  symbols  ( A ),  filled  symbols  (B)y  statistics  is  the  same  as  Fig.  2. 


It  is  important  to  point  out  that,  despite  the  higher  molecular  weight  of  B  (higher 
sedimentation  rate)  in  comparison  to  that  of  A,  the  number  NB  of  B  particles 
becomes  much  larger  than  NA  at  extreme  values  of  H  1.  We  believe  that  this  is 
due  an  implicit  interaction  between  constituents  and  their  source:  the  release  of 
particles  from  the  source  into  the  system  depends  on  their  fraction  in  the  lattice 
causing  a  flow-induced  long-range  hydrodynamic  interaction  among  the  heavier 
constituents.  Both  components  (A,B)  are  driven  upward  by  the  bias  with  equal 
probability.  However,  B  constituents  have  lower  probability  to  leave  the  system  from 
the  top  than  that  of  A  constituents  on  increasing  the  weight  ratio  a  from  1,  which 
results  in  a  higher  concentration  and  therefore  higher  injection  of  B.  One  may  view 
the  molecular  weight  of  particles  as  a  ‘sticky’  long  range  overall  (attractive 

elastic)  interaction  with  the  source  which  increases  with  their  concentration.  This 
results  in  a  non-linear  surge  of  heavier  component  and  decay  of  lighter  component 
accordingly  as  reflected  in  the  flow  response  (to  be  discussed  some  where  else). 

Can  we  quantify  the  dependence  of  the  total  number  (AO  of  these  constituents  as  a 
function  of  hydrostatic  pressure  bias?  In  an  attempt  to  address  this  question,  we  have 
examined  the  variation  of  N  with  the  bias  H.  The  inset  figure  (top)  shows  such  a 
variation  on  a  log-log  scale  which  suggest  a  weak  power-law  decay  N  oc  H~° 14  in  the 
moderate  range  of  H  (see  inset  Fig.  4).  The  plot  on  a  log-normal  scale  does  not  look 
much  different.  The  empirical  dependence  may  be  more  complex  than  just  power-law 
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and  exponential;  perhaps  it  is  a  combination  of  both  with  the  possibility  of  stretched 
exponential  decay  of  N  with  the  bias  H  which  is  difficult  to  identify  with  the  current 
data.  At  large  molecular  weight  ratio  (a  >8-10),  the  non-monotonic  dependence  of 
N  on  H  seem  to  suggest  that  a  non-linear  structure  with  specified  density  can  emerge 
in  such  a  driven  system.  Large-scale  examples  may  include  berm  formation  in 
explosive  detonation  and  in  settling  of  lava  flow  after  a  volcanic  response.  From  the 
transverse  density  profiles  (see  Fig.  3),  one  may  evaluate  the  average  density  dy. 
Variation  of  dy  with  the  molecular  weight  (MB  or  a)  for  different  bias  (see  the  lower 
inset  Fig.  4)  suggests  that  dy  of  A  component  decreases  linearly  while  that  of  the  B 
component  increases  with  a  for  H  >0.2.  At  lower  bias  (H  <  0.1),  transverse  densities 
of  both  components  (A,B)  decays  as  the  sedimentation  dominates. 


4.  Conclusions 

In  summary,  structural  responses  of  a  multi-component  driven  system  are  studied 
by  an  interacting  lattice  gas  model.  Immiscible  components  ( A ,  B)  are  driven  from  a 
source  by  concentration  gradient  and  hydrostatic  pressure  bias  against  gravitational 
sedimentation  in  an  effective  medium.  The  hydrostatic  pressure  bias  H  and  the 
molecular  weight  ratio  are  varied  in  this  open  system.  A  steady-state  structural 
pattern  emerges  at  a  fixed  bias  and  molecular  weight  with  specific  density  profiles  of 
each  component  as  they  continue  their  steady  flow.  In  addition  to  phenomenological 
interactions  among  the  constituents  (Eqs.  (1,2)),  gravity,  and  pressure  bias,  an 
implicit  hydrodynamic  interaction  develops  as  the  injection  of  constituents  in  system 
depends  on  their  lattice  concentration.  The  number  of  each  constituents  reaches  a 
stable  value  (nearly  a  constant)  in  steady  state. 

The  longitudinal  density  profiles  depend  on  bias  and  molecular  weight  with  a 
range  of  patterns  (non-linear,  exponential,  linear)  for  the  dependence  of  density  on 
the  altitude.  The  transverse  density  profiles  show  segregation  with  partial  layering 
for  specific  values  of  bias  and  molecular  weight.  Number  of  B  constituents  (higher 
molecular  weight)  are  higher  than  that  of  A;  the  difference  increases  on  increasing 
the  molecular  weight  of  B.  The  total  number  of  particles  shows  a  weak  dependence 
on  the  bias,  N  cx  H~0A4  for  moderate  range  of  bias  (0.2  <0.8)  for  a  =  1-8.  At 

large  molecular  weight  ratios  (a>  8),  N  becomes  non-monotonic  as  a  function  of  bias 
with  a  minimum  (at  about  FT  >0.83);  the  number  of  A  and  B  constituents  show 
complementary  dependence,  i.e.,  maximum  and  minimum,  respectively,  as  function 
of  bias.  The  large  response  of  the  overall  density  of  heavier  components  at  high  bias 
(H  ->  1)  is  a  signature  of  volcanic  response — a  result  of  the  flow — induced  long- 
range  hydrodynamic  interaction  with  the  source  caused  by  the  concentration 
dependent  in-flux.  Thus,  specific  density  distributions  with  segregation  and  phase 
separation  appear  in  an  open  system — pattern  and  phase  boundaries  can  be 
predicted  at  least  empirically  as  a  function  of  hydrostatic  bias  and  molecular  weight. 
How  realistic  are  these  predictions  remain  to  be  seen  as  the  large-scale  geomarine 
exploration  continue  to  probe  the  morphological  evolution  of  methane  hydrate  and 
related  components  [1-3]. 
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